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REVISITING THE STABILITY OE SPATIALLY HETEROGENEOUS 
PREDATOR-PREY SYSTEMS UNDER EUTROPHICATION 

J. Z. FARKAS, A. YU MOROZOV, E. G. ARASHKEVICH, AND A. NIKISHINA 


Abstract. We employ partial integro-differential equations to model trophic inter¬ 
action in a spatially extended heterogeneous environment. Compared to classical 
reaction-diffusion models, this framework allows us to more realistically describe 
the situation where movement of individuals occurs on a faster time scale than the 
demographic (population) time scale, and we cannot determine population growth 
based on local density. However, most of the results reported so far for such systems 
have only been verified numerically and for a particular choice of model functions, 
which obviously casts doubts about these findings. In this paper, we analyse a class 
of integro-differential predator-prey models with a highly mobile predator in a het¬ 
erogeneous environment, and we reveal the main factors stabilizing such systems. 
In particular, we explore an ecologically relevant case of interactions in a highly 
eutrophic environment, where the prey carrying capacity can be formally set to ’in¬ 
finity’. We investigate two main scenarios: (i) the spatial gradient of the growth rate 
is due to abiotic factors only, and (ii) the local growth rate depends on the global 
density distribution across the environment (e.g. due to non-local self-shading). For 
an arbitrary spatial gradient of the prey growth rate, we analytically investigate 
the possibility of the predator-prey equilibrium in such systems and we explore the 
conditions of stability of this equilibrium. In particular, we demonstrate that for a 
Holling type I (linear) functional response, the predator can stabilize the system at 
low prey density even for an ’unlimited’ carrying capacity. We conclude that the 
interplay between spatial heterogeneity in the prey growth and fast displacement of 
the predator across the habitat works as an efficient stabilizing mechanism. These 
results highlight the generality of the stabilization mechanisms we find in spatially 
structured predator-prey ecological systems in a heterogeneous environment. 


1. Introduction 

Understanding top-down control in food webs has been always a focus of theoreti¬ 
cal ecology, and revealing potential mechanisms which stabilize predator-prey/resource- 
consumer interactions in ecosystems with a high degree of eutrophication (i.e. ecosystems 
that are high in nutrients) is a particularly challenging problem [351 [371IH |55J (H [51 |3T] . 
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Classical theory predicts that a generic predator-prey system will be highly unstable un¬ 
der eutrophic conditions and is expected to exhibit large-amplitude oscillations of species 
densities which will eventually result in population collapse and extinction [121 [H] ■ How¬ 
ever, these theoretical conclusions are often at odds with reality, since there exist a large 
number of examples of ecosystems in which the species densities remain low despite a high 
nutrient supply [ini m [3121 sg. Various mechanisms have been proposed to explain 
such stability of predator-prey interactions in eutrophic environments. In particular, it 
has been suggested that this could be due to interplay between structuring within the 
prey population according to vulnerability or nutrition properties and active food se¬ 
lectivity of the predator (prey switching) |T| [T3 El El Ell IK] . However, stabilization 
mechanisms of this kind cease working for unstructured populations, so a natural ques¬ 
tion is whether efficient top-down control is possible for a population of identical prey 
individuals. 

Spatial heterogeneity of the environment and a non-homogeneous distribution of in¬ 
dividuals across the habitat has been suggested to be an important factor promoting 
top-down control and preventing species extinction [21 El IH H ES] • In recent model 
studies of planktonic systems with a high nutrient supply, it was shown that interplay 
between the gradient of light in the water column and fast vertical displacement of 
herbivorous zooplankton can stabilize the system even for an unlimited nutrient stock 
[301 [27] . The works cited implement a modelling approach which differs from the classical 
reaction-diffusion framework in several important ways. In particular, the new approach 
allows for different time scales for the population growth rate of prey and its predator 
and also takes into account the fact that the rate at which predator disperses across the 
habitat can be substantially faster than its characteristic demographic scale. For exam¬ 
ple, zooplankton herbivores can quickly move vertically along the entire habitat (in the 
euphotic part of water column) a large number of times between reproduction events, 
thus we cannot assign an organism to a particular spatial location as in the case of 
reaction-diffusion models. Thus we should describe the predator population as a whole 
in terms of the total population size, while still taking into account non-homogeneity of 
the distribution of predator individuals across the habitat when considering grazing of 
the prey at each particular location. Ecologically similar scenarios can be found in a 
large number of other ecosystems when modelling, for instance, interactions between ter¬ 
restrial herbivorous mammals and grass, zooplankton and fish, insects and birds, insects 
and their parasitoids, acarine predator-prey systems and even between fish population 
and fishing boats [23 EOl ESI Ell El EH 0 ■ 

Modelling predator-prey interaction in a spatially heterogeneous habitat, where an 
exact location of individuals on a demographic timescale cannot be properly assigned, 
requires the implementation of partial integro-differential equations. However, the com¬ 
plexity of this framework makes it rather difficult to treat the model analytically: all 
previous findings have been obtained by direct numerical simulations of the equations 
for particular parameterisations of the model ingredients [301EZ] • This, obviously, cannot 
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be considered as a rigorous proof of stabilization of the system. Moreover, the results 
can strongly depend on the choice of the per capita population growth rate, which in the 
previous studies was considered to be an exponential function modelling light attenua¬ 
tion with depth [30l[27]. Finally, we should say that the integro-differential equations in 
[80l [27] differ significantly from those used in the modelling of physiologically structured 
populations, see e.g. mill]. For instance, the model we develop here does not contain 
a transport term as standard physiologically structured models do. 

In this paper, we revisit the previous results on stabilization for a class of integro- 
differential predator-prey models, and explore the main factors assuring efficient top- 
down control. We analytically obtain criteria for the existence of a positive stationary 
state assuring the coexistence of both species for an arbitrary gradient of resource dis¬ 
tribution across the environment. Then we investigate the stability of the coexistence 
equilibrium and derive a general characteristic equation, which governs the stability of 
the system. We analytically show that in the case that there is no saturation in the 
functional response (i.e. Holling type I response), the predator is able to stabilize the 
system at low prey density even for an ’unlimited’ carrying capacity. We study how 
the size of the habitat and the spatial gradient of the prey growth affects stabilization 
in the system. We explore two ecologically relevant scenarios: (i) the case when the 
spatial gradient of the growth rate of the prey is due to abiotic factors only, and (ii) 
the case where the growth rate of the prey is affected by the surrounding density of the 
population itself (e.g. effects of self-shading). We find that, interestingly, in case (i) the 
trophic control of the prey is not possible when the population size of the predator varies 
slowly, whereas in case (ii) the predator can efficiently suppress the population growth 
even without changing its own population size. 


2. Modelling framework and biological rationale 


The general predator-prey model with a highly mobile predator is described by the 
following system of partial integro-differential equations m- 




R{h,p{h,t))p{h,t) - f{p{h,t))z{h,t), p{h,0) = po{h), (2.1) 


dZ 

dt 


f{p{h, t))z{h, t) dh - mZ{t), Z{0) = Zq, 


, , k f 

Jo 

Z{t) = ^J z{h,t)dh, P{t) = ^J p{h,t)dh. 


( 2 . 2 ) 


(2.3) 


In the model above p(h,t) and z{h,t) denote the densities of prey and predator at time 
t and at spatial location h, respectively; H is the overall size of the habitat; hence P{t) 
and Z{t) are the average densities of the species at time t; while po and Zq denote the 


initial conditions. The first term in the right hand side of (2.11 describes the dispersal 


of prey across the habitat due to diffusion; R is the function (or functional, in general) 
describing the local per capita growth rate of the prey; / is the local functional response 
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of the predator. The parameters k and m denote the food efficiency of the predator and 
its mortality rate, respectively. In the case where we model plankton interactions, h is 
the vertical coordinate from the water surface and H is the overall depth of the euphotic 
zone. 

Here we impose zero flux boundary conditions = 0 at = 0 and h = H. 

The integro-differential framework (2.1)-(2.3) is quite different from the conventional 
reaction-diffusion modelling setting. Here we consider that the predator movement occurs 
on a much faster timescale than both the prey movement and the population dynamics; 
and assume that the predator can quickly move to any part of the habitat. Thus, we 
cannot assign a particular spatial location to an individual within the predator popula¬ 
tion, since organisms constantly change their locations. The framework still allows us to 
consider the instantaneous spatial distribution of the predator population z{h,t). Note 
that equations (2.1 )-(2.3) do not uniquely specify z{h, t) based on the knowledge of Z{t), 
so we need to impose some further assumptions. A number of empirical observations sug¬ 
gest that the relative distribution of predator between patches often follows the relative 
abundance of food, which is described by 




(2.4) 


The above proportion-based parametrization may not be realistic since it does not take 
into account possible interference between predators |34j . However, in a large number 
of study cases, see e.g. [551 HOI HH HSl ISS], we can still reasonably well approximate 
the instantaneous spatial distribution of predator across the food patches by (2.4). For 
instance, the distribution of herbivorous zooplankton in the water column often follows 
that of phytoplankton, see e.g. EHED]. 

Fig.la shows a typical example of the relative distribution of the biomass of phyto¬ 
plankton (prey) and herbivorous zooplankton (predator) in different layers of the vertical 
water column in the euphotic zone of the sea. Here we split the water column into a num¬ 
ber of layers and plot the values of the spatially averaged biomass of p and z in each layer 
divided by the overall average biomass across the entire column. The data were recently 
collected at 12 stations in the Black Sea for various seasons (2007-2013); further details 
on obtaining the samples are provided in electronic Appendix 1. The figure presents both 
night and day observations. For the night observation (open circles) one can see that the 
data can be roughly described by the straight line with slope approximately equal to 1, 
thus the basic assumption (2.4) holds (standard linear regression gives = 0.85 and 
the slope of the fitting line as 0.96). 

The major difference between the night and the day patterns is that the vertical 
distribution of zooplankton substantially varies due to the phenomenon of regular daily 
vertical migration: during the day time zooplankton is mainly located in deep layers to 
avoid predation by higher visual predators (e.g. fish) and ascend to upper layers at night 
[38] . Note that daytime feeding activity of herbivorous zooplankton is usually low and 
most of the food is consumed during the night time [6]. 
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Fig.lb provides examples of vertical distributions of species and the daily zooplankton 
migration pattern. One can see that the vertical distribution of zooplankton is highly 
variable throughout the 24 hour period. Moreover, zooplankton usually exhibit fast 
short term unsynchronized migrations m. by permanently descending and ascending 
over smaller distances (10-20m). Thus, we cannot assign a particular location to an 
individual zooplankton and we can use the integro-differential framework with an instan¬ 
taneous distribution of the predator. We should note that model (2.1)-(2.3) describes 


the organisms averaged over the day, and therefore excluding daily migration patterns. 
Finally, Fig.lb also provides the profile of chlorophyll a through the column, which is 
a proxy of the distribution of phytoplankton. Its complicated structure is due to high 
heterogeneity in the distribution of the vital resources: light and nutrients (phospho¬ 
rous and nitrogen) through the column, which potentially translates into a complicated 
relation r{h) (cf. [H]). 

We consider that the local functional response is given by the Holling type II function 

[zniiH] 


/(p) = a 


P 


l+pf3’ 


(2.5) 


where the parameters a and /3 denote the attack rate and the inverse saturation food 
density, respectively. In the limiting case when /3 = 0 the functional response becomes 
linear or Holling type I response. Here we intentionally implement a concave downwards 
(i.e. non-sigmoidal) functional response, which is a destabilizing factor in predator-prey 
models in general, see e.g. muni- 

For the sake of simplicity we neglect the effect of diffusion of the prey by setting D = 0. 
In the Discussion section, we shall briefly address the effect that prey diffusion has on 
the results. 

Finally, we need to specify the growth function R. Since we model eutrophic con¬ 
ditions, we consider that this function does not depend on the variation of nutrient 
concentration. We shall explore two main scenarios: (i) the per capita growth rate is a 
function of the abiotic environment only, i.e. it is a local function R = r(h) and (ii) the 
per capita growth rate depends on the population density distribution in some parts of 
the environment R = R{h,p), i.e. it is a non-local term. In the latter case, we have a 
functional, i.e. a function which depends on another function characterizing the overall 
spatial distribution of p. Note that in both cases, we neglect natural mortality of the 
prey compared to the predation rate. 

In case (i) the equations read (Model 1) 


dp 

m 


{h,t) 


= r{h)p{h,t) 


dt 


(t) = k 


m 

HP{tf 


" pHh.t) 


0 1 + I3pih,t) 


dh — mZ{t), Z{0) = Zq. 


( 2 . 6 ) 

(2.7) 
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For the moment we do not specify the function r. We only assume that it is a smooth 
and non-negative function. 

In case (ii), we shall investigate the particular scenario of self-shading in the popula¬ 
tion. In this complex scenario, the supply of a resource such as light decreases in space 
via a function of the local population density. In planktonic communities this may be 
shading of light for phytoplankton residing in deeper layers by the organisms dwelling in 
upper layers. For a non-planktonic system, this may be for instance a gradual depletion 
of a nutrient as it flows through part of the environment occupied by the prey and is 
gradually consumed by the organisms. In this case, the local rate of decrease of resource 
density can be assumed to be a function of the density of consumers. It is natural to 
assume that the available resource quantity is a decreasing function of the total num¬ 
ber of consumers p(x,t) dx that are above/upstream an individual at depth/point h. 
In particular, we assume that the quantity of the available resource is an exponentially 
decreasing function of the number of consumers that contribute to the shading. 

The equations read (Model 2) 


dp 


(/i, t) = To exp u J p{x,t) p{h,t)— 


^(t) = k^ar 

dt^’ HP{t) Jo l+(3p{h,t) 


p'^ih,t) 


( 2 . 8 ) 

(2.9) 


dh — mZ(t), Z{0) = Zq, 


where u is the coefficient of self-shading. In both of the models above the average densities 


P{t) and Z{t) are given by equation (2.3). 

In this paper, we are mostly interested in the existence and stability of the positive 
stationary state of (2.6)-(2.7) and (2.8|- (|2.9[ ). It is worth mentioning that a spatially 
homogeneous predator-prey model (2.1 )-( [2l^ with perfect mixing (i.e. for r{h) = const 
and z{h) = const) will be globally unstable [4], thus the main question is whether it 
is possible to stabilize this system by introducing spatial heterogeneity (i.e. structu¬ 
ring the population) and allowing a non-homogeneous distribution of predators. It is 
also worth noting that the term describing prey growth in equation (2.8) has infinite 


dimensional nonlinearity, and that structured population models with such nonlinearity 
are notoriously difficult to analyse, see e.g. 0, US- 


3. Results 

3.1. Model 1: The local growth rate of the prey is a function of abiotic factors 
only. 


3.1.1. Existence of a non-trivial (coexistence) steady state. 
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We look for a strictly positive steady state of model (2.6)-(2.7). We assume 

that r is a strictly positive continuous function with 

max W(h)} = tm- 
he[0.H] 

We characterise the existence of the positive steady state in the following proposition. 

Proposition 3.1. Assume that 


ak > to/3 


holds. Then model (2.61-(2.71 has a unique positive steady state. 
Proof. The steady state equations read 

\ P^ijk) L ^ In 


Zf, — 

p* = 


P* 1 + /3p:t{h) ' 

Pl{h) 


i-H 


ak 

mH Jq l + ^p*(/i) 


d/i, 


rH 


Pt{h) dh. 


For positive values {Pt,Z^.) which satisfy 

aZ* > PtmT*, 


we obtain from equation (3.111 


p,(h) = ^ 


"{h) 


ajp - r{h)l3 


, hG[0,H]. 


(3.10) 

(3.11) 

(3.12) 

(3.13) 

(3.14) 

(3.15) 


We then substitute the steady state p* given by formula (3.15) into equations (3.12)- 
(3.13). After simplihcation we obtain: 


n 1 jr 

H Jq ajr - r(h)/3 

(3.16) 

- ^ r d/I 

(3.17) 

mH Jo of. - r(/i)/3 ' 


From equation (3.17) we obtain: 

r I - r{h)l3 


—- 


mH 


0 


+ 4 




k 

/3mH , 

from which we obtain 


fiajr-r{h)(3 Pajr-r{h)l3) 

k aZ 

TO 13 


'(/i) 


dh 


^\^h)dh- 


Zs, — 


k 


1 


rH 


ak — m/3 H Jq 


r{h) dh. 


(3.18) 

(3.19) 
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We substitute this expression into equation (3.16) to obtain 

f 

10 


1 = 








ak 

r{h) Ah— - - — r{h)l3P^ 


Ah =; F{P,). 


(3.20) 


/o 


ak — mf3 


Also note that (3.14) together with (3.19) require the existence of a (unique) solution of 
(3.20) in the interval 


ak 


1 1 


ak — m/3 /Atm H 


J r{h) d/i^ . 


Note that F{0) = < 1, and that 


(3.21) 


F/Psf) — y + 00 , 


as 




ak 


1 1 


cH 


ak — m/A /Atm H Jq 


r{h) Ah 


(3.22) 


Hence there is a solution of equation (3.20) in the required interval (3.21). It is also 


shown that on the interval (3.21) the function F is strictly monotonically increasing, 
hence the solution is unique, and the proof is complete. □ 

Note that, we may relax the assumption of strict positivity of r, and allow it to vanish 


at some values of h. Then equation (3.15) shows that the steady state p* will also vanish 


at those values of h, but Proposition 3.1 still holds. 


3.1.2. Stability of the coexistence steady state. 

The stability of the coexistence equilibrium is determined by the roots of a character¬ 
istic equation, which we derive in the current subsection. The linearisation around the 
coexistence steady state (p*,Z*) of model (2.6)-(2.7) is computed as 

dAp 


dt 


-{h,t) =r{h)Ap{h,t) 


Pli.h) 


P*(l -I- /Ap^/h)) 


AZit) - 


pl{h)Z^ 


AAZ , , . s ak . , 

-(t) = - mAZ{t) + -—AZ{t) 


PH1 + Pp*ih)) 
Pl{h) 


^P(t) + Ap(ft, t) 


cH 


di 


HP. 


1 -\- Pp.{h) 




P*(l -I- l3pt{h)Y 

(3.23) 

pH 


pi{h) 


+ 


akZ^ p*(h)(2 -I- j3p^,{h)) 


lo A+Pp*ih) 

cH 


Ah 


HP. 


{1 + /3p.{h)y 


Ap{h,t)Ah, AP/t) = ^ f Ap(h,t)Ah. 

H Jo 

(3.24) 
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We look for the solution in the form of Ap{h, t) = w{h) exp(At) and AZ{t) = U exp(At). 
This leads to the eigenvalue problem {W is the spatial average of w) 

p*(/i)(2 + Pp*{h))' 


A — r{h) + ( 


P* (1 +/3p*(/i))2 

Pl{h) 


P*(l + Pp*{h)) 


U- 


w{h) = 

pUh) 


XU = -mU 
ak 

+ -T7- 


Pl{h) 


p2 1 + I3p^{h) 
pl{h)Z. 


W, 


u- 


H Jo \P*il +I3p*ih)) P2(l 


(3.25) 




Using (3.15) we recast equation (3.25) as follows 


X + r{h) ( 1- ^-r{h) 
Z* a 


j{h) = 


P*(l + l3p^(h)Y 

P^r'^{h) 

aZ^-P^Pr{h)" Z^{aZ^ - P^I3r{h)) 


Jh) 


-W- 


(3.26) 

U. (3.27) 


Note that (3.14) implies that aZ* — P*r(/i) > 0, for/i € [0,P]. Next we multiply equation 
(3.25) by then integrate it from 0 to P and add it to equation (3.26) to obtain: 


XkW + XU = kW- mU, 


where 




pH 


We re-arrange equation (|3.27 1 as 


;(h) = 


Hh) 


-W- 


r{h)w{h) d/i. 


\h) 


(3.28) 

(3.29) 


P, 


fci(A, h) {aZ^ — P^j5r{h)) ki{X, h) {aZ^ — P^pr{h)) Z^ 


U, hG[0,H], 


(3.30) 


where 


ki{X,h) = X + r{h) ( 1 - ^-r{h) 

z* a 


h e [ 0 ,P]. 


Note that A:i(A, •) 7 ^ 0 for A S C, Re(A) > 0. 


Next we integrate equation (3.30) from 0 to H and divide by H to obtain 


W = k2{X)W - k2{\)^U, 


(3.31) 


where 


1 


r^{h) 


Ah. 


Iq fci(A, h) {aZ^. — P^j3r{h)) 

We also multiply equation (3.30) by and integrate from 0 to P to obtain 

_ , , , P* 
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where 


1 /■" 


r^{h) 


0 fci(A, h) {aZ^ - P^(3r(h)) 


Ah. 


Next we substitute into equation (3.28) the expression for W given by (3.32). This 


together with equation (3.31), yields a homogeneous system of two equations for the 
variables {W, U) as follows. 


0 =k (ksiX) -X)W+ (^-k^ksiX) -X-m]U, 
0=(l-fe(A))lT + fc2(A)^t/. 


(3.33) 


For any value of A S C with Re(A) > 0 a non-trivial solution of system (3.33) yields a 
non-trivial solution w via equation (3.30). On the other hand, any eigenvalue A with 


Re(A) > 0 necessarily satisfies equations (3.33) for same non-trivial {W, U). 
We summarize our result in the following theorem: 


Theorem 3.2. The non-trivial steady state (p*, Z^) of model (2.6)-(2.7) is locally asymp¬ 


totically stable if every solution X of the characteristic equation below has negative real 
part. 

K{X) := k^ksiX) - k 2 {X) (^A (^1 + k^^ +m^ +X + m = 0, (3.34) 

where k 2 and k^ are defined earlier. On the other hand, the steady state is unstable if 


there is at least one solution of (3.34) with positive real part. 


In general, the derived characteristic equation (3.34) is complicated and determining 


its roots might require extensive numerical simulations which might take as much time 
as directly solving the system of differential equations. However, this equation can be 
rather useful for construction of a Hopf bifurcation curve: one can set A = (6 is a 

positive real number) and solve the resulting equation. In the Discussion section we shall 
provide an example of construction of such a curve for a particular parametrization of 
the function r. 

Using the characteristic equation above we are able to address the stability of the 
system in the case when the demographic timescale of the predator is much slower than 
that of the prey. In this case, we can assume the population size of the predator to be 
constant. The result is given by the following remark. 


Remark 3.3 The positive steady state (p*, Z^) is always unstable under perturbations 
in the first component (in the prey population) only. We consider perturbations in the 
p-subspace only. Perturbations in the first component only lead to the characteristic 
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equation 


1 f r^{h) 

fc 2 (A) ki{X,h){aZ^ - Ptj5r{h))^^ 

_ 1 r^{h) 

H Jo (^x + rih) (l - ^ £r(/i)) ) {aZ, - pP.r{h)) 


dh. 


(3.35) 


It is shown that fc 2 ( 0 ) > 1, and clearly ^2 is a monotone decreasing function of A for A > 0. 
Also we have that lim fc 2 (A) = 0, hence there exists a unique positive eigenvalue, and 

A—)-+oo 

the steady state is unstable. 


Remark 3.4 Another important conclusion we may draw from (3.34) is that for a 
homogeneous environment (i.e. r{h) = const) the system is always unstable provided 
/3 > 0, since all eigenvalues will have positive real parts (we do not show here the 
corresponding characteristic equation for the sake of simplicity). 


3.1.3. Stability in the special case /3 = 0. 

Next we investigate an important special case of (3.34) where we have a Holling 
type I functional response, i.e. there is no saturation rate in predation, ,5 = 0. The 
characteristic equation reduces to 


P h 

K{X) =m + X+^- / 


r^{h) 


0 {XPr{h))aZ, 


dh 


- A 1 + fc 


Zi, 


i-H 


\h) 


1 


H Jo {X + r{h))aZg 


■ d/i = 0 , 


(3.36) 


which can be simplified by recalling the expressions for the stationary densities of species 

(m-|-A)A r{h)dh 


(3.16), (3.17) and (3.19). We obtain 

r'^{h)dh 


K{X) =m — 


2mX 


J^r^{h)dhJo A-br(/i) J^r{h)dhJo A + r(/i) 


= 0 , 


(3.37) 


Let us first consider the case of a small environment, i.e. small H. We use the Taylor 
expansion of (3.37). We assume that the derivative of r{h) at h = 0 does not vanish: 
r{h) « ro + Rh, where i? 7 ^ 0. The characteristic equation becomes 


rom + A^ 1 TT 1 


-R^X 


Tom — 3Xm — SrpA + A^ 


+ 0{H^). 


A + ro ’ 2''''‘"(A + ro)2*‘ 12ro^'' " (A + ro)^ 

(3.38) 

This expression can be re-written in the polynomial form (A —rp) by dropping the 

term 0{H^) 


0 — uq -h aiX -h a2X^ -h a^X^ -h n 4 A^, (3.39) 

where oq = rgm, oi = m(2rQ — i7^i?^/12 -|- HroR/2), 02 = mro + Tq — (rg — m)RH/2 — 
R^{3m + 3rg)RV(12rg), 03 = 2rg - RH/2 - R^H^/{12ro), 04 = 1. 
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We use the Routh - Hurwitz stability criterion for the coefhcients of the above poly¬ 
nomial: an > 0, 0203 > 04 O 1 , 01 O 2 O 3 > 040 ^ -I- aoo|. One can see that for a sufficiently 
small H: o„ > 0 (n = 1,2,3,4), moreover 


01O2O3 — 0403 — ooOg = mrl{ro + m)R^H '^/2 + 0 {H^) > 0, 
0203 — 0401 = 2 rQ -I- 0{H) > 0. 


(3.40) 

(3.41) 


Thus, for small H the polynomial characteristic equation (3.39) has eigenvalues with 


only negative real parts and, consequently, the initial characteristic equation (3.37) also 


exhibits stability for small H. Note that stability is sustained in the case when we 
consider a parabolic approximation r{h) « tq -f Rh -\- Rih? , where i? 7^ 0,i?i 7^ 0 (the 
need for considering such an approximation comes from the fact that in the second order 
expansion ( 3.38[ ) the quadratic term i?i will be present. We should stress that here we are 
allowed to drop 0{H^) and consider (3.39) since all the roots of (3.39) are perturbations 


of the solutions with H = 0 such that the maximal order of perturbation with respect to 
H does not exceed 2, i.e. it is less than 3 (we do not show the corresponding expressions 
X{H) of solutions of (3.39) for the sake of brevity). 


Now let us suggest that for some critical H > 0 the stability of (3.37) changes and the 


real part of an eigenvalue becomes positive. This can occur at a Hopf bifurcation point 
and the condition of stability change is A = z&, where 6 is a positive real number. We 
substitute X = ib into the complex characteristic equation and obtain two equations for 
its real and imaginary parts. After some re-arrangement the equation for b becomes (see 
electronic Appendix 2 for the exact derivation) 


- 2 




1 + 


,2 r(h)dh 

^ Jo b^+r^(h) 
r{h)dh 


pH r(h)dh 

Jo b^+r^(h) 
/(f r{h)dh 


= 0 . 


(3.42) 


We substitute a particular parametrization r{h) to verify whether or not (3.42) has a real 


value solution. In the case (3.42) does not have real solutions, the steady state of (|2.6|)- 


(2.7) will be always locally stable for /3 = 0 because for small H the eigenvalues always 


have non-zero negative real parts and for any H > 0 they will not vanish (otherwise we 


would have a solution with A = ib, in which case b would solve (3.42)). 


However, it is still an open question whether or not the equilibrium will remain sta¬ 
ble for (3 = 0 for any arbitrary positive function r{h). In our investigation, we checked 
several concrete families of ecologically relevant and mathematically tractable functions 
including linear, parabolic, cubic, exponential and sinusoidal functions given, respec¬ 
tively, by r{h) = qq + aih, rih) = oq -I- 02 / 1 ^, r{h) = oq -f 03 / 1 ^, r{h) = aoexp(a 4 /i) 
and r{h) = ao(l -I- esm{ujh)), where ai,u},e are parameters. For those functions we can 
obtain analytical expressions for the corresponding integral and numerically check the 


possibility to find a solution of (3.42), by varying the parameters. For all these functions 


we find that (3.42) has no real solution, thus the equilibrium is always stable (we do not 


show the corresponding cumbersome expressions and graphs here for the sake of brevity). 
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3.2. Model 2: local growth rate of prey is affected by self-shading. 

3.2.1. Existence of a non-trivial (coexistence) steady state. 

We look for the coexistence equilibrium (p*, Z*) of model (2.8l-(2.9). The steady state 
equations read 


roexp vp{h,t)dt?j = 


P*{h) 


+ j3p*{h)' 


rH 


m = 


P. = 


Pl{h) 


ak 

P*H Jq l + /3p*(/i) 
1 


d/i, 


f-H 


H , 

From the first equation of the system we obtain 

rh 


In(ro) — u / p{h,t)dh = \n 

Jo 


We differentiate the above expression to obtain 

1 


p*(/i) dh. 

Z^, p*{h) 

P* 1 + l3p*{h) 


Integration of (3.47) gives 


-vh + C =- 


Pl{h){l + I3p*{h)) 

1 + 


1 

p*{h) 


■pin 


P*{h) 


(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 


where C is a constant of integration. The constant C can be found by setting h = 0 and 
using the expression for p* (0) from ( |3.43 1 

P*ro 


p*(0) = 


(3.49) 


aZ^ — cqP^P 

After substituting C and some simplification, we obtain the expression for the stationary 
state Pt{h): 


—vh + 


r[)PifP — aZ* 

roP* 


■pin 


( aZ^, 


\roP* 


1 


P*{h) 


■ P In 


1 + Pp*{h) 
p*{h) 


(3.50) 


The above expression is implicit since it is impossible to solve it analytically for pif(h) 
except in the particular case /3 = 0; however, from (3.47) one can conclude that p*(/i) 
is a decreasing function of h. The maximal density through the column is attained at 
h = 0 and given by (3.49); the minimal density is achieved a.t h = H. From (3.43) we 

P*ro exp{—HiyP^,) 


obtain the expression for p*(i7) 

P*{H) = 


aZ^ - tq exp(-iJuP*)P*/3’ 


(3.51) 
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We use (3.51) and (3.50) to get the following relation between P, and 

0 = — exp{PtHv)aZt + pP^Hv + HvtqP^ + aZ^. 


(3.52) 


Next we consider the stationary equation (3.44) for the predator and replace the 


integration variable h with p{h) using (3.47). 

Pl{h) 


ak 

Kh 


rH 


1 + I3p*{h) 


dh = — 


ak 


P.{H) 


P^Hv (1+ /3p*(/i))' 


■ dp = 


ak 


P^Hvp 


p,(H) 


l + /?Pjp.(0) 

(3.53) 


We are allowed to make the above change of variables since p{h) is a monotonically 


decreasing function. We substitute the values ofp*(0) andp*(i7) from (3.49) and (3.51) 


respectively, to obtain after simplification another identity linking P* and Z* 

(1 — exp(—P*Pu))roA: 


Zs, — 


mHv 


(3.54) 


Finally, we substitute the expression for Z^ from (3.54) into (3.52) to obtain the equation 
for P* which reads 


f{P,) = l3P^Hu + P,Hiyro + 


akrn 


{exp{—P:tHiy) — l)(exp(P*Pi/) — 1) = 0, (3.55) 


mHv 

This is a transcendental equation; however, we can make exhaustive conclusions about 
the existence of a positive solution of (3.55). At low P*, the function /(P*) is linear 


and always has a positive slope. At large P*, /(P*) becomes negative (the exponential 
term becomes dominant), thus there should be at least one point where this function 


vanishes. Additionally, equation (3.55) should be considered together with the condition 


p*(h) > 0, which is equivalent to p*(0) > 0 since p*(h) is a decreasing function. From 
(3.49) and (3.54) it follows that the condition p*(0) > 0 is equivalent to 

fi{P*) = P^PmHv + ak exp{—P^Hiy) — ak > 0, (3.56) 


It is easy to see that (3.56) will be satisfied for some P* > 0 if and only if ak — /3m > 0. 


Further, by computing the second derivative of f{P*) one can show that for ak — 13m > 0 
it is always negative, thus there can be only one positive solution of (3.55). Hence, 


combining the above properties of /(P*) and /i(P*), we find that the existence of the 


nontrivial steady state of the system (2.8)-(2.9) is given by the following 


Proposition 3.5. Model (2.8)-(2.9) admits a unique positive steady state (p*,Z*) in the 
case ak — /3m > 0 and the unique positive solution of (3.55), which would always exist 


under the mentioned condition, satisfies the inequality (3.56). 


It is easy to see that for /3 I the conditions of the proposition are satisfied, thus there 
always exists a non-trivial steady state. Let us consider a gradual increase in /3 while still 
keeping ak — /3m > 0. The root of /(P*) = 0 will increase, whereas the root of /i(P*) = 0 
will decrease. At a certain /3* both roots coincide, which is an unavoidable event since 
the root of /i(P*) = 0 will approach zero when ak = /3m. Thus for all /3 < (3*, the system 
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has a non-trivial coexistence equilibrium whereas for /3 > /?*, a positive equilibrium is 
impossible. The exact value of /?* can be found from the condition /(P*) = /i(P*) = 0 
which results in a cumbersome analytical expression (not shown here). 


3.2.2. Stability of the non-trivial (coexistence) steady state of Model 2. 

The linearisation around the positive steady state of (2.8)-(2.9) results in 

a rather cumbersome characteristic equation. In this paper we focus on an important 
particular case, where there is no saturation in the functional response of predator, i.e. 
/? = 0. The linearized system becomes 


(h, t) =ro exp y vp^(h)dh^ Ap(/i, t) — J Ap{x,t)dx 


dAp 


dAZ afcZ* 

dt ^ ~ HP? ^ 


f 2p^{h)Ap{h,t)dh. 

JO JQ 


(3.57) 

(3.58) 


We look for the solution in the form of Ap{h, f) = w{h) exp(At) and AZit) = U exp(At). 
We substitute the expressions for the stationary solutions (3.48) (/3 = 0) and ( |3.54 ) to 
arrive at the following eigenvalue problem (W is the spatial average of w) 


}(h)X = 


Cro 
hu + C 

Cro 


m- 


hv + C Jq 


}{x)da 




U\ = 


kCro 

H 




cH 


1 


oH 


'o (hu + cr 


dh + 2 


}{h) 


hv + C 


dh 


(3.59) 

(3.60) 


We re-arrange the first equation of the above system 

ch 


w{h)\{hv-\-CY = — CrQ{hvC)w{h) — vCtq f w{x)dx — Cro f ^ , (3.61) 

Jo J 


We differentiate both sides of (3.61) and obtain 
' dw 


dh 


dw , 


[hi/ -I- C) -I- 2w{h)v X{hi^ C) =Cro —vw{h) —-\- C) — vw{x) (3.62) 


dh 


After simplifying the above equation we get 

((hv + C)X + vC) (hv + C)'^=- 2vw(h) ((hv + C)X + vC) 

dh 


(3.63) 
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We suggest that (hu + C)X + vC ^ 0 and derive the following equation for the eigen¬ 
function w{h) 


[hv + ~ 2i/w{h) 

Integration of this linear equation gives the following eigenfunction 

= (h/+cr ^ 

Next, we substitute this function into system 
X = -ro - Cro [ — -— 


C/A = 


kCro 

H 


H 


-a"'. 


1 


cH 


’o [hiy + cy 


where the value of W is obtained from 

rH 


I rH I 

^ J H j 


Ah+ 2 


dh 


1 


'o {hu + Cr 


dh 


We simplify equation (3.67) 


UX = 


kro 


{hv + CYC 


{hv + CY C{Hv + C) ■ 


- +2C + HV 

P* 


(3.64) 

(3.65) 

(3.66) 

(3.67) 

(3.68) 

(3.69) 


finally, we substitute the expressions for the stationary value P* into (3.66) and (3.69) 
to obtain the characteristic equation 

0 = A^-|-(l — Qi)roA -|- Q 2 , 

where 

Hv 


Q2 — 


Qi — 


krn 


{C + Hv) ln{{C + Hv)/C) 


< 1, 


(3.70) 

(3.71) 


Hv 


• + 2C + Hv] > 0. 


{hv + CYCZ^ V ln{{C + Hv)/C) ■ — ■ “-y — (3-72) 

It is easy to see that characteristic equation ( |3.70 ) always has negative real parts, thus, 
the linearized system (3.57)-(3.58) is always stable. 

We should also say that unlike the model without self-saturation the number of eigen¬ 
values is hnite and there is a single eigenfunction. This can be explained by the fact that 
separation of variables provides us with only the solution for large values of time and 
does not describe transient regimes at small and intermediate times. 

The above result can be summarized as the following 


Theorem 3.6. The steady state {pf,ZY of model ( |2.8[ )-(2.9) is always locally asymptot¬ 
ically stable for j3 = 0 (Holling type I functional response). 
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Remark 3.7 One can prove that in a similar way the stability of equation (3.57) 
governing the dynamics of prey for a constant (equilibrium) density of predator, i.e. for 
Z{t) = const or U{t) = 0. The eigenvalue of the system is given by 


A = - (1 - Qi)ro < 0. 


(3.73) 


This signifies that a spatial regrouping of the predator population (without variation 
in the total population size) alone can control the prey population in the case of self¬ 
shading. 


4. Discussion and summary of results 


In this paper, we revisit predator-prey interaction in ecosystems with a high nutrient 
load and explore the role that spatial heterogeneity of the environment and high predator 
mobility have in stability of such ecosystems. Our analytical investigations of the model 
predict that interplay between these two factors results in stabilization of an otherwise 
globally unstable predator-prey system. An important assumption for the realization of 
the given stabilization scenario is that the predator preferentially feeds at locations with 


high prey density (see relationship (2.4)). Since spatial gradients in a species’ growth rate 


across their habitats is a typical ecological situation, and many predators/consumers are 
highly mobile and feed mainly in patches of high food density [291 EDI EH EH 0 [35] (see 
also figure 1 of the current paper), we claim here that the reported mechanism should 
be rather typical for marine, freshwater and terrestrial ecosystems under eutrophication. 
Our findings also provide a new and robust mechanism to resolve the famous ’paradox 
of enrichment’, which is a long standing enigma in theoretical ecology [HiiiiiiiiiTiEa 
[331143]. 

Our models are formulated as integro-differential equations, which allows for the incor¬ 
poration of different time scales: fast spatial movement of predator and slow population 
dynamics of both species. This makes the current framework more effective than clas¬ 
sical reaction-diffusion models for modelling trophic interaction in small-sized habitats, 
since individuals cannot be assigned to a particular spatial location in this case while 
reaction-diffusion equations contain only local terms. 

Here we explore a generic predator-prey model with an infinitely large carrying ca¬ 
pacity of prey (i.e. no resource limitation) and consider two main scenarios: (i) the prey 
growth rate gradient is fixed (Model 1) and the local growth rate of prey is affected by 
self-shading (Model 2). The main mathematical outcomes of our study are the following: 

(i) For each model (Model 1 and Model 2) we find the condition for the existence of the 
positive stationary state through which prey and predator can coexist for an arbitrary 
gradient of the growth rate r{h) > 0. Interestingly, including self-shading usually restricts 
the conditions of coexistence (cf. Propositions 3.1 and 3.5). 

(ii) For Model 1 we have derived the general characteristic equation (3.34) governing 


the stability of the coexistence equilibrium. Using this equation, one can construct the 
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Hopf bifurcation curve that separates stable and unstable dynamics without numerically 
integrating the underlying partial integro-differential equations (see an example below). 

(hi) Analysis of the characteristic equation (3.34) shows that for a homogeneous en¬ 
vironment, in which r(h) = const, the system is always unstable provided there is sat¬ 
uration in the functional response (Remark 3.4). Thus, spatial heterogeneity of the 
environment is a necessary condition for stability of the coexistence equilibrium. 

(iv) For both Model 1 and Model 2, in the case of a linear (Holling type I) functional 
response of the predator, we analytically proved that the coexistence equilibrium is always 
stable. Since the system is structurally stable, the stability of the equilibrium is also 
guaranteed for sufficiently small /3 > 0, i.e. for a Holling type II functional response with 
low saturation. 

(v) Comparison of Model 1 and Model 2 shows that stabilization in the system with 
dynamical self-shading (Model 2) is potentially a stronger mechanism than the scenario, 
where the spatial gradient in the prey growth rate is fixed (Model 1). For instance, 
without saturation in the functional response (/3 = 0) the stabilization does not require 
changing the predator biomass: variation in the spatial distribution of a fixed number of 
predators can suppress the prey outbreak and re-establish the equilibrium (cf. Remarks 
3.3 and 3.7). 

In the case of a fixed prey growth rate gradient, the general characteristic equation 
(3.34) allows us to easily compute the Hopf bifurcation curve to direct numerical sim¬ 
ulation of the underlaying integro-differential equations (2.6)-(2.7). As an illustrative 
example, here we have constructed the Hopf bifurcation curve for the system with a 
linear variation of the growth rate r{h) = oq + aih (see electronic Appendix 3 for detail). 
Fig.2 shows a set of bifurcation curves in the diagram oi and (3 constructed for a gradu¬ 
ally increasing size H of the habitat. In each case, the stability region is located below 
the curve. One can see that increased the saturation in food consumption, signified by 
a large /3, usually acts as a destabilizing factor, but the influence of the gradient of the 
growth rate ai is more complicated: small and large ai result in destabilization of the 
system whereas stabilization occurs for some intermediate oi. An increase in the size of 
the habitat can also either destabilize or stabilize the system. This result is somewhat 
counterintuitive since in earlier works it was shown that for an exponential parametriza- 
tion of r{h), an increase in the spatial gradient of the growth function or the size of the 
habitat will always facilitate stabilization of the system EHIIIT]. 

Our analytical investigation confirms the numerical results reported earlier in [Ml HZ] 
about the possibility of the stabilization of systems with ’unlimited carrying capacity’, 
and extend the previous findings to more general parametrizations of the spatial gradient 
of r{h). In particular, we proved that for a small habitat size H, any arbitrary non- 
homogeneous dependence r(/i) will stabilize the predator-prey interactions (see Section 
3.1.3). On the other hand, our work shows that using some other parametrizations, for 
example a linear function, can provide some qualitatively different prediction compared 
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to the exponential dependence considered in [33 HZ]. This underlines the danger of 
sticking to particular parametrizations when modelling predator-prey systems. 

Throughout the paper we have neglected the diffusion of the prey population by set¬ 
ting Z? = 0 in equation (2.1). Including diffusion D > 0 will make the model more 
realistic, but the equations become analytically untractable. Although one can still per¬ 
form linearization techniques, it becomes impossible to obtain the resultant characteristic 
equation due to the presence of the second order derivative with respect to x, thus the 
only way of exploring the system is by direct numerical simulation. For small values of D 
the main conclusions of our paper obtained for D — 0 will hold due to the structural sta¬ 
bility of the model equations. Interestingly, for some ecosystems (e.g. planktonic systems, 
[33 [31]) the diffusion coefficient D > 0 is relatively small (compared to the interaction 
terms) so the results obtained should be close to those with D = 0 . However, for a very 
large diffusion coefficient corresponding to a highly mobile prey, the system becomes ho¬ 
mogeneous in space and is reduced to the classical Rosenzweig-MacArthur model, which 
is always globally unstable for any 13 > 0. Thus, the system will be destabilized for a 
certain critical D, which can be found numerically. 

Finally, we should mentioned the open questions that our study raises. For instance, it 
would be an important extension to determine the stability or otherwise of the coexistence 
equilibrium of Model I in the case r{h) is an arbitrary function for a Holling type I 
functional response (see equation (3.42)). Based on preliminary results obtained for a 
small habitat size H, we hypothesize that for ,5 = 0 the Model I should be stable for an 
arbitrary r(h), but it remains to prove this conjecture rigorously. Another challenging 
issue is exploring the system stability for a more complicated relationship between z(/i, t) 
and p{h,t), by taking into account possible interference between predators |34j . For 
instance, this can be represented by setting z{h) = Z, /i > 0 [33] ■ Another interesting 
extension would be to investigate analytically the model combining self-shading of the 
population with a fixed gradient of local growth rate and combining Model I and Model 
2. From the biological point of view it might even be natural to consider situations when 
the growth function R (or r) can take negative values, i.e. the population of prey actually 
declines locally. In this more general case, model (2.6)-(2.7l will still be well-posed, and 
solutions starting in the positive cone would remain non-negative for all times. At the 
same time equation (3.11) shows that there will be no strictly positive steady states of 
the model and this might influence the previous stability results. Models with locally 
negative growth rate of prey could therefore possess quite rich dynamics, and such models 
should definitely be a matter of further investigation. 
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B 


mg C/m3 mg Chl/m3 rng C/m3 mg Chl/m3 



Figure 1 . (A) Relative distribution of phytoplankton and zooplankton 
biomass across the water column collected at 12 stations in the Black 
Sea. Details on plankton sampling are provided in electronic Appendix 
1. For each station, the whole euphotic zone is split into a number of 
layers. Figure shows the relative phytoplankton biomass plotted against 
the relative zooplankton biomass in the same layer. Here Pi and Zi are 
the average densities of species in layer i; P and Z are the overall aver¬ 
age densities. 

(B) Typical example of vertical distribution of zooplankton [mgC/m^) 
and phytoplankton (mgChl/rrt') during night (bottom left panel) and 
daytime (bttom right panel) collected at station #5. The profile of 
phytoplankton distribution is averaged over 24 h. High variation in 
zooplankton distribution through the day is due to the diel vertical mi¬ 
gration (for more explanation see text). 
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Figure 2. Hopf bifurcation curves constructed for model (2.6)-(2.7) for 
the case the growth rate is a linear function of space r(/i) = uq + aih. 

The curves are plotted using the derived characteristic equation (3.34), 
for details see Appendix 3. Various curves are obtained for different 
size H of the habitat and each of them separates the stability region 
(below the curve) from the instability region (above the curve). The 
other model parameters are m = 0.1; oq = 1; fc = 1 ; = 2 . 
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